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ABSTRACT 

By fitting a flexible stellar anisotropy model to the observed surface brightness and line-of-sight velocity 
dispersion profiles of Draco we derive a sequence of cosmologically plausible two-component (stars + dark 
matter) models for this galaxy. The models are consistent with all the available observations and can have 
either cuspy Navarro-Frenk- White or flat-cored dark matter density profiles. The dark matter halos either 
formed relatively recently (at z ~ 2 ... 7) and are massive (up to 5 x 10'^ Mq), or formed before the end of the 
reionization of the universe (z ~ 7 ... 11) and are less massive (down to ^ 7 x 10^ Mq). Our results thus support 
either of the two popular solutions of the "missing satellites" problem of A cold dark matter cosmology - that 
dwarf spheroidals are either very massive, or very old. We carry out high-resolution simulations of the tidal 
evolution of our two-component Draco models in the potential of the Milky Way. The results of our simulations 
suggest that the observable properties of Draco have not been appreciably affected by the Galactic tides after 
10 Gyr of evolution. We rule out Draco being a "tidal dwarf" - a tidally disrupted dwarf galaxy. Almost radial 
Draco orbits (with the pericentric distance < 15 kpc) are also ruled out by our analysis. The case of a harmonic 
dark matter core can be consistent with observations only for a very limited choice of Draco orbits (with the 
apocentric-to-pericentric distances ratio of < 2.5). 

Subject headings: galaxies: individual (Draco dwarf spheroidal) — galaxies: evolution — methods: A^-body 
simulations — stellar dynamics 



1. INTRODUCTION 

Galactic dwarf spheroidal galaxies (dSphs) are intrigu- 
ing objects with a deceptively simple appearance - roughly 
spheroidal shape, no gas, and no recent star formation. Due 
to their closeness, these galaxies are studied in significantly 
more detail than other external galaxies (see the review of 
lMateJl998 ). Despite that, the nature of dSphs and their place 
in the larger cosmological picture is not well understood. 

The first wave of enhanced interest in dSphs took off af- 
ter the pioneer work of A aronson (J983). Based only on the 
three stars in Draco with measured line-of-sight velocities, 
he boldly claimed that Draco can be significantly dark mat- 
ter (DM) dominated. The fact that the stellar velocity disper- 
sion in dSphs is significantly larger than what would follow 
from the virial theorem for the luminous mass was later con- 
firmed wit h much larger studies. F or example, the latest com- 
pilation of Wilkinson et alJ J2004h contained 207 Draco stars 
with measured Une-of-sight velocities. Similar results were 
obtained for other dSphs as well - both for the Milky Way 
and M31 satellites. There were attempts to explain the large 
stellar velocity dispersion in Galact ic dSphs withou t invoking 
the DM hypothesis. Most notablv. FKroupal J 1 997h presented 
a model where dSphs are considered to be "tidal dwarfs" - 
virtually unbound stellar streams from dwarf galaxies tidally 
disrupted in the Milky Way potential. The model of Kroupa 
(0,997) appears to be not applicable to Draco as it is unable 
to reproduce the narrow observed horizontal branch width of 
this dwarf (Klessen et al. 2003). In this paper we will present 
additional evidence against Draco being a tidal dwarf. The 
overall consensus now appears to be that dSphs do contain 
significant amounts of DM and are hence the smallest known 
objects with (indirectly) detected DM. This makes them very 
interesting objects, as they can be an important test bench for 
modern cosmological models and for the theories of DM. 

More recently, the interest in dSphs was rejuvenated after 
the realization that simple (DM-only) cosmological models 



overpredict the number of the Milky Way satellite ga laxies 
by 1- 2 orders of magnitude (.Klvpin et al, 1999: Moore'et alJ 
Il999h . This was coined the "missing satellite problem" of 
cosmology. The original analysis was done under the assump- 
tion that DM in dSphs has the same spatial extent as stars. Re- 
laxation of the above assumption led to a suggestion that only 
the most massive subhalos predicted to populate a Galaxy- 
sized halo mana ged to form stars, w ith the rest of the subha- 
los staying dark dStoehr et al.l2002h . Another way of solving 
the missing satellites problem is to assume that only the oldest 
subhalos formed stars, with the star formation in the younger 
subhalos being suppressed by the metagalactic ionizing radi- 
ation after the r eionization of the un iverse was accomplished 
around z ~ 6.5 (iBullock et al."2000'). In reality, both mecha- 
nisms could have realized (Ricotti & Gnedin 2005). 

To discriminate between the two above solutions of the 
missing satellites problem and to place Galactic dSphs in 
the right cosmological context one has to know the global 
parameters of these dwarfs, and most importantly, their to- 
tal DM extent and mass. The traditional approach is to as- 
sume that the dwarf is in equilibrium (thus ignoring the pos- 
sible impact of the Galactic tides), and to solve the Jeans 
equation to infer the density profile of the DM halo based 
on the observed surface brightness and line-of-sight veloc- 
ity dispersion profiles. As the proper motions of individual 
stars in dSphs are not known, one has to resort to making 
certain assumptions about the anisotropy in the stellar veloc- 
ities. Due to a well known degeneracy between the assumed 
stellar anisotropy and inferred total enclosed mass there are 
many solutions to the Jeans equation which are consistent 
with the observations. Another limitation of the above ap- 
proach is that DM can be traced only within the stellar body 
extent of the dwarf, so no conclusion can be made about the 
total mass of the galaxy. It is also not clear at what dis- 
tance from the dwarf's center the virial equilibrium assump- 
tion breaks down due to Galactic tides. Some work has been 
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done on the impa ct of tidal forces on t he structure of Galac- 
tic satellites Te.g. iQh. Lin & Aarsethlfl 995: Piatek & PrvoS 
1995HHavashi et alJl2003HKazantzidis. Maeorrian. & Moorg 



20041) . where it was clearly demonstrated that the tidal strip- 
ping and shocking is a complex dynamic process. The defi- 
ciency of the above work is in using single-component mod- 
els for the dwarf galaxies, which made it impossible to di- 
rectly compare the results of the numerical simulations with 
the observations. In general, stars in dSphs are distributed 
differently from DM, so they also behave differently in reac- 
tion to the external tides. To correctly describe the observable 
manifestations of the Galactic tides in dSphs it is essential to 
use two-component (stars + DM) dwarf models jRead et al.l 
|2005). 

In this paper we place joint constraints on the global prop- 
erties of Draco (one of the best studied dSphs) by (1) using 
cosmological predictions for the properties of DM halos, (2) 
developing very flexible stellar anisotropy model for dSphs, 
and (3) running a large set of high-resolution simulations of 
the evolution of two-component dwarfs in the Galactic tidal 
field. We derive a sequence of cosmologically plausible mod- 
els for Draco which are consistent with all the observed struc- 
tural and kinematical properties of this dwarf. Our results are 
consistent with either of the two above solutions for the miss- 
ing satellite problem. 

2. GLOBAL CONSTRAINTS ON DRACO DM HALO 
PROPERTIES 

We consider two types of DM halo s: theoretically moti- 
vated EivOToI^e^I&mi^ j 1 9971 hereafter NFW) halos 
with a 7 = — 1 density cu sp at the center, and observationally 
motivated lBurkertI i 1 99.51) halos, which have a flat core: 
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(l + r/r,)[l + (rA,)2] 
Here and r, are the scaling density and radius. At large 
distances from the center both halos have the same asymptotic 
density slope of 7 = -3. 

Analysis of cosmological A^-body simulations showed that 
the concentration c = rvir/^j of low-mass DM halos (with the 
virial mass rnvk- = 10*^ . . . lO" M0) has a log-normal distribu- 
tion with the mean 
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and dis persion 0.14 dex fe ullock et al. 2001, with the correc- 
tion of Sternberg. McKe e. & Wolfirell2002l) . Here ryir is the 
vi rial radius of the halo and z is the redshift. 

Sternberg et al. (2002) showed that the four dwarf galax- 
ies with a Burkert DM halo profile studied by Burkert (.1995.) 
have the same dependence of the DM halo scaling density ps 
on the scaling radius r, as do the NFW halos in cosmological 
simulations. This result was obtained for z = 0. We assume 
that it holds true for other redshifts as well, and use the equa- 
tion (0 to find concentrations of both NF W and B urkert halos. 

Using the formula of Sheth & Torme one can es- 

timate the comoving number density of DM halos per unit 
Inmvir and per standard deviation in concentration: 
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Here Vc is the number of standard deviations from the mean 
concentration, v = (0.707/5)'''^(5(z), where S(z) is the criti- 
cal overdensity for spherical collapse extrapolated linearly to 
z = 0, p„, is the present day mass density of the universe, and 
S is the variance of the primordial density field on mass scale 
OTvii- extrapolated linearly to z = 0. To estimate the above pa- 
rameters we follow Mashchenko, Couchman, & Sills (2005). 
Throughout this paper we assume a flat ACDM cosmology 
and use the following values of the cosmological parameters: 
n„ = 0.27, nb = 0.044, // = 71 km s"' Mpc"', and erg = 0.84 
(Spergeletal. 2003). 

It is interesting that one can derive quite general and (stel- 
lar) model-independent constraints on the properties of the 
Draco DM halo by combining the available observational data 
on this galaxy with the predictions of cosmology. Throughout 
this paper we assume that the distance to Draco is 82 ± 6 kpc 
(Mateo 1998). It is convenient to consider different DM halo 
constraints in the plane of two scaling parameters, ps and r^. 
We summarize the global constraints in Figure^ 

The most obvious constraint is that the Draco DM halo 
should have formed before the bulk of the Draco stars formed. 
We assume that the first star burst in Draco took place at least 
10 Gyr ago ( Mateo- 1998.) . or at z ^ 1 .8. In Figure[T] the areas 
to the right of the dashed Unes marked "z = 1.8" correspond 
to halos older than 10 Gyr. 

The next constraint, F ^ Fmin, comes from the requirement 
for DM halos to be abundant enough to explain the observed 
number (^ 20) of dwarf spheroidal galax ies in the Local 
Group. Following Mashchenko et aQ (120051) . we adopt F^m = 
0.01 Mpc"^. As the function F does not explicitly depend 
on Ps and (it depends on mvi, and z), in Mashchenko et aQ 
(2005) we designed a numerical method for finding the most 
probable combination of (mvir,z) corresponding to given val- 
ues of ips,f's)- As a by-product we also obtain the correspond- 
ing value of i/f. The areas below the solid lines in Figure ^ 
correspond to DM halos which were abundant enough to have 
been progenitors of dwarf spheroidals in the Local Group. 

Assuming that the tidal field of the Milky Way has not 
perturbed significantly the stars in the outskirts of Draco 
out to a distance of rout ~ 1-2 kpc from its center (the last 
observed point in the D raco surface brightness profile of 
'Odenkirche n et al]l200lL see Figure |2^), the third constraint 
can be written as rvir ^ rout- (It is hard to imagine stars form- 
ing beyond the virial radius of its DM halo). The halos in the 
areas above the dash-dotted lines in Figure^satisfy the above 
criterion. 

In Figurel^we plot the observed line-of-sight velocity dis- 
persion profile for Draco (from Wilkinson et al. 2004). It has 
been noted that in Draco and Ursa Minor the galacti c out - 
skirts appear to be kinematically cold ( Wilkinson et al.l2004^ . 
There appears to be no good explanation for this phenomena. 
Given the fact that in the case of Draco the only (outermost) 
radial bin which is "cold" contains only 6 stars with mea- 
sured line-of-sight velocities, and is marginally (at ^ 1 sigma 
level) consistent with the global average for Draco, (crios) = 
9.5 km s"', we decided to exclude the last bin from our anal- 
ysis. Our assumption is that at the distance of ri = 0.74 kpc 
from the Draco center (half-way between the two last points 
in FigureO the Une-of-sight velocity dispersion for Draco is 
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Fig. 2. — (a) Stellai' surface density profiles for Draco. Tlie observed profile of 'Odenkirchen et alj 120011 their sample S2) is shown as solid circles with 
one-sigma error bars. The best fitting profiles for Plummer-like model with a = 7, Plummer model, and King model are shown as solid, dotted, and dashed lines, 
respectively, (b) Stellar anisotropy profile Ti(r) for the best-fitting Draco model with NFW DM halo (a = 7, logp, = 7.2, logr, = 0.45, A = 1, rjo = 0.9, rji = —0.7) 
is shown as a solid line. For comparison, anisotropy profiles are shown for other values of A: 0.5 (dotted line), 3 (short-dashed line), and 5 (long-dashed line). 



roughly equal to the global average': aiosiri) ^ 9.5 km s"' 
(the circle with a cross in Figurel^. At this distance the stel- 
lar density in Draco is steeply declining (see Figure [2^). As 
a result, most of stars observed at the projected distance ri 
from the Draco center are located roughly in the plane of the 
sky at the spatial distance ri from the center of the dwarf. 
We can then use crios('"i) as a lower limit of Draco circular 
velocity at this distance, Vd- Indeed, if one considers the 
extreme of purely radial orbits, the line-of-sight velocity dis- 
persion at large distance from the galactic center will be close 
to zero. In the opposite extreme of purely circular orbits, one 



' This is consistent with the results of lMunoz et alj 120051) who observed 
Draco to have crjos 10 km s"' out to a distance of ~ 1 kpc from the 



can show that a\os becomes approximately equal to the circu- 
lar velocity at this distance. We plot the solution of the equa- 
tion Vc.i = 9.5 km s~' as long-dashed lines in Figure [T] The 
areas above these lines correspond to halos which satisfy the 
above criterion. The result we derived here is approximate, 
but of model-independent nature. We present more accurate 
(but also more model-dependent) treatment in § |3l where we 
fit stellar models to all the reliable observed ctios points in Fig- 
ure|3] 

The last constraint is that the virial mass of the Draco 
halo is somewhere between "ridiculously" low and high val- 
ues 10^ and 10'' M0 (the area between two dotted Unes in 
Figure The lower limit is even lower, by a factor 2 - 3 
than the classical "mass follows light" estimates ( lMateoll998t 
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Fig. 3. — Lin e-of-sight stellar velocit y dispersion profile for Draco. Obser- 
vational data of lWilkinson et al] llOO-ft are shown as circles with one-sigma 
eiTor bars. The last point which we consider to be unreliable is shown as 
an empty circle with error bars. The horizontal dotted line marks the global 
CTlos = 9.5 km s"' for Draco. The (Tjos profiles for our 8 best-fitting stellar 
models are shown as solid lines. The big circle with a cross shows our as- 
sumption for crios at the distance of ri = 0.74 kpc from the Draco center. 



lOdenkirchen et alJl2001h . The upper limit corresponds to a 
satellite which would very quickly spiral in to the center of 
the Milky Way due to dynamical friction. As one can see, 
the last constraint does not add any new information to our 
exclusion plots in Figure 

From comparison of Figures^ and[n) one can see that all 
the constraints except for the fourth one (Vc i ^ 9.5 km s"') are 
very similar for both NFW and Burkert DM density profiles. 
This is intimately linked to our assumption that for a given 
virial mass niyir and redshift z both types of halos have the 
same concentration - given by equation (|3}- This automati- 
cally makes the scaling radii r, equal in both cases. At the 
same time, one can show that the scaling densities for NFW 
and Burkert halos are related through 



ln(l+c)-c/(l+c) 



p5, Burkert _ 2 

Pi.NFW ln( 1 + c) + [ln( 1 + c^)] /2 - arctanc ' 



(5) 



which is close to unity for realistic halos (the ratio changes 
from 0.966 . . .0.921 for c = 3.5 . . . 10). The fourth constraint, 
unlike the rest of the criteria, does not deal with a global prop- 
erty of a halo (such as myir and z or their derivatives - rvir and 
F). Instead, it deals with the average DM density within a 
certain fixed radius - which can be dramatically different for 
cuspy NFW and flat-cored Burkert models. 

The shaded areas in Figure[Ocorrespond to DM halos which 
satisfy all of the above global constraints. As one can see, the 
Draco halo parameters are not particularly well constrained, 
especially in the case of the NFW profile. Still, one can 
make a few interesting observations. Firstly, the most restric- 
tive (and useful) constraints are F ^ 0.01 Mpc~^ and Vc.i ^ 
9.5 km s"'. Secondly, there is plenty of room for Draco to 
have formed before the end of the reionization of the universe 
at z ~ 6.5 (shaded areas to the right of the dashed lines marked 
"z = 6.5" in Figure^. Thirdly, we can derive the range of pos- 
sible values of different halo parameters. For NFW halos, our 



exclusion plot implies that lognjvii = 7.7 . . . 10.7, z = 1 .8 ... 11, 
= 0. 16 ... 7.8 kpc, and log pj = 7.0 . . . 9. 1 . (Here units for 
mvh and ps are Mq and M0 kpc""*, respectively.) For Burkert 
halos, logmvi, = 7.8 . . . 10.5, z = 3.9 ... 1 1, = 0.17 .. .5.0 kpc, 
and log Ps = 7 A. . .9.1. Interestingly, for both NFW and Burk- 
ert cases, Draco could not have formed before z^ 1 1, or more 
than 13.2 Gyrs ago. The reason for that is that at larger red- 
shifts DM halos with the virial radii ^1.2 kpc become too 
rare to correspond to a typical dwarf spheroidal galaxy in the 
Local Group. 

The main purpose of generating the exclusion plots in Fig- 
ure n was to significantly reduce the computational burden 
in the next step of our analysis, described in the next sec- 
tion, where we use the whole observed line-of-sight velocity 
dispersion profile along with the surface brightness profile to 
find the best-fitting stellar models and to further reduce the 
uncertainty in {ps,f's) values for the Draco DM halo. 

3. STELLAR MODEL 

The equilibrium state of a spherically symmetric stellar sys- 
tem can described by the Jeans equation. 



1 d(p^a^) ^ 2 
dr r 



dr 



(6) 



iBinnev & Tremain3ll987ll) . which is obtained by taking the 
first velocity moment of the collisionless Boltzmann equation. 
Here r is the distance from the center of the system, is 
the stellar density, cr, and cr, are the one-dimensional stellar 
radial and tangential velocity dispersions, respectively, and $ 
is the total gravitational potential (due to stars and DM). The 
radial gradient of the gravitational potential can be calculated 
as d^/dr = G[m(r) + m^(r)]/ p-, where m{r) and m*(r) are the 
enclosed DM and stellar masses, respectively, and G is the 
gravitational constant. From equations (QJI^Jl we derived 

m(r) = 47rr^p,[ln(l+jt:)-jt:/(l+jc)] (NFW), (7) 

m(r) = 27rrj Pi [ln( 1 + x) + ( 1 /2) ln( 1 + jc^) - arctanjc] 

(Burkert). (8) 

Here x = r/rs. We neglect impact of baryons on DM distri- 
bution, as in Draco the stellar density is more than an order 
of magnitude lower than DM density even at the center of the 
galaxy. 

The Plummer density profile. 



p,=po[l + (r/bf] 



-5/2 



(9) 



jBinnev & Tremain^ll987^ . is used sometimes to describe 
simple spherically symmetric stellar systems, such as glob- 
ular clusters and dwarf spheroidal galaxies. It has a core of 
size b and a power-law envelope with the slope 7 = -5. We 
found that a "generalized Plummer law". 



p. = p„[l + (r//7)y"/', 
which has a surface density profile of the form 



(10) 



(11) 



provides much be t ter fit to the Draco star count profile of 
lOdenkirchen et all i2001l) than the Plummer model and the 
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theoretical iKind lfT966ll model, if one chooses a = 7 (see Fig- 
ure|2t) ^. Here a is an integer number ^ 2 and R is the pro- 
jected distance from the center of the system. For a ^ 4, the 
total stellar mass can be calculated as 



Mh.= 



a-3 

For a = 7, the stellar enclosed mass is 



mjr) ■■ 



Airpor^ 5 + 2{r/bf 
15 (l + rV/72)5/2' 



(12) 



(13) 



and po = 15So/( 16fe). The fitting of eq uation ( fTH to the 
Draco profile of lOdenkirchen et"aD (1200 A gave po = 108 x 
10^ M© kpc"^ and b = 0.349 kpc. (We assumed that the stellar 
V-band mass-to-light ratio of Draco stars is T = 1.32, which 
is an average value of the Salpeter and composite model esti- 
mates of Mateo et al. 1998.) 

Traditionally, in equation (|6} one uses a,- and (3=1- (jf/cr^ 
instead of cr,. and cr,. The anisotropy parameter f3 is equal 
to -oo, 0, and 1 for purely tangential (circular), isotropic, and 
purely radial stellar orbits. We advocate a different anisotropy 
parameter. 



(14) 



Unlike (3, parameter ry is symmetric: it is equal to -1, 0, and 
1 for circular, isotropic, and purely radial orbits. The equa- 
tions connecting f3 and 77 are 77 = (3/(2- (5) and (3 = 2r?/(l +77). 
Another useful expression is = al{l-rj) / (l + ri). The Jeans 
equation (|6} can now be rewritten as 



1 dip^af.) ^ 4?7 < 



dr 



(15) 



p* dr 1 + ry r 

As there are two unknown functions in equations (|6} or 
M5\ . (T,(r) and cr;(r) (or (3[r\, or r][r\), one customarily as- 
sumes the shape of the anisotropy profile and then solves the 
ordinary differential equation for <Jr{r). The boundary condi- 
tion is ar{r — > 00) = 0. Traditional choices for the anisotropy 
profile are (1) (3 = constant (with the special case of /3 = 0, 
or isotropic stellar orbits), and (2) the Osipkov-Merritt profile 
rOsiDkov,1979:Men-itt1985l). 



/?= l + {ralry 



(16) 



In the latter case, the stellar system is isotropic at the cen- 
ter, reaches (3 = 0.5 at r = r^, and becomes purely radially 
anisotropic in infinity. Unfortunately, the two above choices 
are very limited. We propose instead a much more flexible 
anisotropy profile. 



ry = 770 + (771-770) l-(p*//Oo) 



(17) 



which is applicable to systems with a flat core (such as gen- 
eralized Plummer model or King model). Here A is a positive 
number of order of unity, and 7/0 and -qi are asymptotic values 
of the anisotropy parameter -q for r and r ^ 00, respec- 
tively. As one can see, the profile in equation Ml\ does not 



^ For the Draco star count of lWilkinson et all 120041) the best-fitting value 
of a is 6. 



explicitly depend on r, like the Osipkov-Merritt profile. In- 
stead, it depends on the stellar density p,. We believe it is a 
reasonable approach, as in many realistic stellar systems the 
same dynamical processes shape simultaneously both density 
and anisotropy profiles. The examples are the collapse of ini- 
tially homogeneous warm stellar sphere (va n Albada 1982t 
Mashchenko & Sills 2005a) and the expansion of a newly 
formed stellar system in a spherical galaxy with a DM halo 
after removal of the leftov er gas by the feedback mechanisms 
('Mashchenko et al. 2005). In both cases, the stellar orbits be- 
come increasingly radially anisotropic in the outskirts of the 
relaxed system, where the density is steeply declining. 

Both (3 = constant and the Osipkov-Merritt profiles can be 
considered as special cases of our more general expression in 
equation (I17> . Indeed, fixing 7/0 = 771 would correspond to the 
case of (3 = constant; using the generalized Plummer density 
profile from equation (I10> and setting A = a/2 produces the 
Osipkov-Merritt profile with rg = bj \fl. 

The parameter A in equation (I17t controls how sensitive the 
anisotropy parameter is to changes in density. To illustrate 
this effect, in Figure|2j3 we show anisotropy profiles for Draco 
with 770 = 0.9, 771 = -0.7, and A = 0.5, 1,3,5. One can see that 
by varying A by a factor of a few a large range of possible 
anisotropy profiles is produced. 

We designed a numerical algorithm to find values of the 
parameters controlling the shape of the anisotropy function 
f] ivo, Vi' ™d ^) which would minimize the deviation of 
the simulated line-of-sight velocity dispersion profile from the 
observed one (circles with error bars in Figure |3l — we ex- 
cluded from our analysis the last point as an unreliable one) 
for given ps, r,, and the halo profile (NFW or Burkert). The 
procedure consists of the following six steps. 

(1) We choose values of logps and logr, from a grid with 
spacings of 0.45 dex and 0.15 dex, respectively. (The values 
of the increments were chosen to lead to a factor of two in- 
crease in the halo virial mass.) The reference point of the grid 
is log ps = 9 and log rj = 0. Only those grid points are consid- 
ered which are lying within the shaded zones in Figure ^ In 
the Burkert case, we also considered a few points lying out- 
side of the shaded area in attempt to bracket the point with 
the best x^- All these grid points are shown as circles (ei- 
ther empty or filled) in Figure |3 Overall, we considered 35 
different DM halo models. 

(2) For each of the DM halo models, we consider 
1764 different combinations of anisotropy shape parameters 
r;o = -1,-0.9,-0.8,...,!, 77, =-1,-0.9,-0.8,...,!, and A = 
0.5, 1,3,5. For each combination, we solve the Jeans equa- 
tion ( I15t numerically with the anisotropy, stellar density, en- 
closed DM mass, and enclosed stellar mass profiles given by 
equations illi . (|9jl, (EHD^ ™d il3i . respectively. We adopt 
a = 7, po = 1 -08 X 10^ M0 kpc-\ and b = 0.349 kpc. The so- 
lution of the Jeans equation is the radial velocity dispersion 
profile (Tr{f)- 

(3) For each of the above ^ 60,000 models we generate a 
spherically symmetric A^-body stellar model, with the number 
of particles = 10,000. Stellar particles are distributed ran- 
domly, with the density profile given by equation (|9jl. Each 
stellar particle is assigned random values of the radial and 
two tangential components of the velocity vector: Vr, Vg, and 
Vtj,. All three components are assumed to have a Gaussian 
distribution, with the dispersions o'r(r), a,{r), and (J,{r), re- 
spectively. This is an approximate method of generating a 
close-to-equilibrium A^-body system with an arbitrary density 
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Fig. 4. — Best-fitting stellar models for different Draco DM halo parameters p, and r,. Solid (empty) circles correspond to models which have < 9.5 
(X^ > 9.5) between the modeled and observed line-of-sigh velocity dispersion profiles. Thick solid lines mark the areas where all the global constraints on 
(Ps,r,) are satisfied (see FigureQ. Thin solid lines correspond to isothermal stellar models with a = 9.5 km s"'. Short-dashed lines correspond to halos formed 
at z = 6.5. Long-dashed l ines con'espond to halos with the analytical tidal radius equal to 0.85 kpc for the extreme values of the orbital pericentric distance, 
Rp = 2.5 and 70 kpc (see S I4.1I . Circle with a cross marks the best-fitting Draco model from lMashchenko et al J i2005 ) for the Draco surface brightness profile of 
lOdenkirchen et al. |2001). 



and anisotropy profiles. The accurate method would involve 
numerically calculating the distribution function, which is a 
very computationally expensive procedure. This would ren- 
der our approach unfeasible. 

(4) We calculate the line-of-sight velocity dispersion pro- 
file for the generated A^-body stellar models integrated over 
the same projected radial bins as in the observed pro- 
file of Wilkinson et al. (2004) (the edges of their bins are 
0,0.12,0.24,0.36,0.48,0.60,0.72 kpc - Jan Kleyna, private 
communication). For this we use the projection method of 
[Mashchenko & Sills (2005a, Appendix B) where we explic- 
itly use the spherical symmetry of our stellar system. 

(5) We calculate difference between the observed line- 
of-sight velocity dispersion profile (the six reliable points in 
Figure |3} and the modeled one. As the observational error 
bars are asymmetric, for calculating we use the appropriate 
one-sided value of the standard deviation depending on which 
side of the observed point the model point is located. 

(6) For each DM halo model, we choose one of 1764 mod- 
els, with different r/o, rji , and A, which produces the lowest 
value of x^. 

For most DM halo models, we could not find a good fit 
to the observed line-of-sight velocity dispersion profile (with 
some models having > 100): all modeled crios('") points 
were either well above or well below the observed ones. Only 
a few models produced x^ < 9.5 (solid circles in Figure]^. It 
is interesting that for both NFW and Burkert cases the best- 
fitting models follow very closely a sequence of isothermal 
stellar models with the total one-dimensional stellar velocity 
dispersion dtot = [(cr,^+2(T,^)/3]'''^ = 9.5 km s~' (thin solid lines 
in Figure 0}. We obtained the isothermal solutions by solving 
the appropriate Jeans equation, 

;^^^+7^"'-''-)=-i^' ^^^^ 



with the boundary condition (7^(0) = (Ttot- For isothermal sys- 
tems the usual Jeans equation requirement CT,.(r — oo) = 
does not hold in general case, and at some radius r^ax the so- 
lution breaks down when aj becomes negative. The isother- 
mal model lines shown in Figure |4] were obtained by find- 
ing the value of which would maximize r^ax for given p^. 
Typically r^^x ^ 1 kpc, and only for the NFW model with 
Ps = 10^ Mq kpc"^ did it drop down to 0.6 kpc. 

The closeness of our best-fitting models to isothermal mod- 
els does not imply that the isothermal models are the best 
ones. We calculated x^ differences between the observed and 
modeled line-of-sight velocity dispersion profiles for a few 
isothermal models, and they were substantially worse than for 
our best-fitting models. We also plotted atoiir) profiles for the 
best-fitting models, and they were not isothermal. The ex- 
planation for the closeness of our best-fitting models to the 
isothermal ones is in the virial theorem. Given that the line- 
of-sight velocity dispersion profile and the surface brightness 
profile are fixed, the virial theorem implies that all realistic 
models should stay close to a certain line in the (p.s, Ts) plane. 

From Figure |4] one can see that the best-fitting models are 
well bracketed inside the zones were all the global constrains 
on Ps and from §|2are satisfied. In other words, we have a 
good agreement between the global constraints and the de- 
tailed a\os{r) analysis constraints, which were derived in a 
very different fashion. The exception is the Burkert halos with 
log p5 = 8 . 1 , where x^ is slightly improving when moving up- 
wards into the area with F < 0.01 Mpc"^. This is explained 
by the fact that at this value of log the isothermal sequence 
is making a sharp upward turn. 

We list the parameters for the best-fitting models in Tabled 
We show only the models with x^ < 9.5 located within the 
globally constrained zones. As for the Burkert halos with 
log Ps = 8 .55 we have two comparably good models, we chose 
the one which is closer to the isothermal sequence. 

Analysis of Tabled shows that a comparably good crio.s(?") 
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TABLE 1 
Best-fitting stellar models 



Input Parameters 


Derived Parameters 


Model 


Halo 


log At 


logr, 


A'dm 


EDM 


A 


m m X" "ivir Z C ''vir 




logF 






Mq kpc"' 


kpc 




pc 




Mq kpc 




Mpc"' 


Nl 


NFW 


7.20 


0.45 


10*" 


60 


1 0.9 -0.7 5.4 4.6 x lO"* 2.46 5.55 15.6 


-0.68 


-1.09 


N2 


NFW 


7.65 


0.00 


3.16 X 10^ 


28 


1 


0.7 -0.8 6.6 5.2 x 10* 4.50 4.80 4.80 


-0.23 


-0.15 


N3 


NFW 


8.10 


-0.30 


105 


36 


1 


0.8 -0.9 6.1 1.8 X 10** 7.10 4.55 2.28 


0.54 


0.08 


N4 


NFW 


8.55 


-0.60 


105 


20 


1 


0.6 -1.0 8.9 6.9x 10' 9.45 5.14 1.29 


1.47 


-0.11 


N5 


NFW 


9.00 


-0.75 


10= 


19 


1 


1.0 -1.0 8.6 8.4 X 10' 10.7 6.92 1.23 


2.80 


-1.62 


Bl 


Burkert 


8.10 


0.15 


10« 


37 


1 


0.7 0.0 6.6 4.5 X lO' 6.74 4.98 7.03 


1.48 


-1.93 


B2 


Burkert 


8.55 


-0.45 


3.16 X 10' 


30 


1 


0.7 -0.6 6.5 2.1 X 10* 9.64 5.17 1.84 


1.82 


-0.97 


B3 


Burkert 


9.00 


-0.75 


105 


19 


1 


0.8 -0.9 6.2 9.1 X 10' 11.0 6.91 1.23 


2.90 


-1.83 



Note. — Here A'dm and eqm are the number of the DM particles and the gravitational softening length for DM particles in A'-body simulations (see § 14.21 . 



fit can be obtained for a very large range in virial masses, 
from ~ 7 X 10^ to - 5 X 10'' Mq, for both NFW and Burk- 
ert halos. Despite very large difference in DM masses and 
density profiles, all the best-fitting stellar models have com- 
parable values of the anisotropy parameters: A = 1 (in other 
words, anisotropy 77 is a linear function of the stellar density), 
r/o ~ 0.8, and r/i ^ ... - 1 . 

In Figure |3l we show the line-of-sight velocity dispersion 
profiles for the 8 best-fitting models from Tabled integrated 
over the same projected radial bins as the observed profile. All 
the models correctly reproduce the observational trend of a 
slight increase in (j\os with radius. However, two models with 
the largest virial mass (both NFW and Burkert) produce pro- 
files which are rising rather steeply at the last measure d point. 
This c ould present a problem if the result of Mun oz et al.l 
that the ctios profile in the outskirts of Draco is almost 
flat, is confirmed with a larger sample of stars. For less mas- 
sive models, the line-of-sight velocity dispersion profiles are 
levelling off at the last measured point, which is more in line 
with the results of Munoz et al. (2005). 

Assuming that the best-fitting stellar models follow closely 
the isothermal sequence, from Figure |3 we can derive new 
(slightly better than in the previous section) constraints on 
Draco's DM halo parameters: logmyir — 7.9 . . .9.7 (for both 
NFW and Burkert halos), z = 1 .8 ... 10, = 0.21 ... 3. 1 kpc, 
logp, = 7.0. . .8.8 (for NFW halos), and z = 6.8 ... 11, = 
0.18 .. . 1.4 kpc, \ogps = 8.1 . . .9.0 (for Burkert halos). The 
interesting result is that if cosmological halos have a Burkert- 
like flat-cored DM density proflles, then Draco should have 
formed before the end of the reionization of the universe at 
z-6.5. 

4. EVOLUTION IN THE MILKY WAY POTENTIAL 

4. 1 . Possible Orbits in the Galactic Potential 

To carry out tidal stripping simulations for Galactic satel- 
lites, it is of principal importance to know reasonably well the 
shape of the gravitational potential of the Milky Way. One 
popular Milky Way model often used to calculate orbits of 
Galactic globu lar clusters and dwarf spheroidals is that of 
Llohnston et alJ ( il999 ). This model consists of three compo- 
nents: (1) a disk represented by Mivamoto & Nagai ( 1975) 
potential with the mass 1.0 x 10" M0, radial scale-length 
6 .5 kpc, and scale- height 0.26 kpc, (2) a spherical bulge with 
a|Hernauist ( 1990) potential, mass 3.4 x 10'° M©, and scale- 
length 0.7 kpc, and (3) an isothermal halo with cr = 128 km s"' 
and a flat core of size 12.0 kpc. The model was designed to 



reproduce the observed flat Galactic rotation curve between 1 
and 30 kpc. 

There are two major disadvantages of the above model for 
our work. Firstly, it is axisymmetric, which adds an additional 
degree of freedom to our already multi-dimensional problem. 
(Orbits in axisymmetric potentials depend on all three compo- 
nents of the space velocity vector of the satellite, whereas in 
spherical potentials orbits depend only on the radial, V,., and 
tangential, V,, sp ace velocity components.) Secondly, the DM 
halo of Johnston et alJ J 19991) is an isothermal sphere with the 
rotation curve asymptotically approaching Vc= 181 km s"' as 
r 00, whereas in the cosmological halos the rotation curves 
are declining in the outskirts. 

We decided to use a simple static NFW potential, 

^ = -At:GqMIMI+RIR,)IR, (19) 

as the Milky Way model for our project. The scaling ra- 
dius Rs and density Qs of an NFW halo can be determined 
from the virial mass Myir and concentration C. These quan- 
tities are still not very well known for the Milky Way. Tra- 
ditionally, one uses different Galactic objects (stars, gas, 
globular clusters, dwarfs spheroidals) with known line-of- 
sight velocity and, in some cases, proper motion as kine- 
matical tracers of the Galactic potential, with different au- 
thors obtaining quite different results. The favored model of 
^Ivpin. Zhao. & Somerville (2002), for example, has a virial 
mass of 10'^ M0 and concentration 10... 17. Other recent 
NFW halo ba^ed models have Mvh- = (0.7-1.7) x 10'^ Mq, 
C = 5 ... 12 fCard one & Sereno 2005), andMvir = (0.6-2.0) x 
10^2 Mq, C= 18 (Battaglia et aL .2005) . Non-NFW mod- 
els of [Sakamoto et al. (2003) give larger values for the to- 
tal mass of the Milky Way: Mvir = (1.5-3.0) x 10'^ Mq 
if Leo I is gravitationally bound to the Galaxy, and Mvi, = 
(1.1-2.2) X lO'^M© if not. 

We adopted an intermediate value for the Galactic virial 
mass, Mvir = 1.5 x 10'^ M©. The median concentration of 
cosmological halos with such mass at z = is C = 13.2 
( Bullock et al. 2001). The halo scaling parameters are then 
Rs = 22.6 kpc and = 6.0 x lO*" Mq kpc"^ and the virial 
radius is /?vir = 298 kpc. 

The radial tidal acceleration -d^^/dr^ of o ur NFW halo is 
comparable to that of the composite model of lJohnston et alJ 
( 1999). As one can see in Figure|5j the tidal acceleration pro- 
file for the NFW halo is located between the two extreme pro- 
files for the Johnston et al. C1999 ) model (in the Galactic plane 
and along the polar axis) down to a radius of 10 kpc. 
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TABLE 2 
Draco's orbits 



Fig. 5. — Radial tidal acceleration profiles for the Milky Wa y models. Solid 
line correspond to our NFW model. The composite model of ljohnston et alj 
IT99 9) is shown as dashed (in the Galactic plane) and dotted (in the polar 
direction) lines. 



We obtained the pericentric and apocentric distances, Rp 
and Ra, for closed Draco's orbits in the potential given by 
equation ( |19> by solving numerically the following non-linear 
equation (£innev & Tremaine 1987): 



1 2{^{R)- 



$(7?o)]-y/-y,2 



2t/2 



= 0, 



(20) 



where Rq, V,-, and V, are the Draco's current distance from 
the Galactic center, radial velocity, and tangential velocity, 
respectively. If proper motion is known, the two space ve- 
locity vector components Vr and V, can be calculated using 
the procedure outlined in Appendix. The radial orbital pe- 
riod is obtained by solving numerically the following integral 
JBinnev & Tremainell987h : 



P = 2 



dR 



(21) 



IScholz & IrwinI Jl994h published the only available mea- 
surement of the proper motion of Draco: jiaCosS = 0.6 ± 
0.5 mas yr"' and jig = 1-1 ±0.5 mas yr"'. (We adopted the 
larger value for the uncertainty from the text of their paper.) 
When we started this project, we hoped that the measure- 
ments of Scholz & Irwin ( 1994) could be used to place use- 
ful constraints on possible Draco's orbits in the Milky Way 
potential. Unfortunately, this is not the case. This can be 
seen from Figure |6 k, which shows the proper motion mea- 
surements of iScholz & IrwinI J 19941) and the locus of the pos- 
sible proper motion vectors resulting in a closed orbit in 
the Milky Way potential given by equation M9\ . (We as- 
sumed that the halo density drops to zero beyond the virial 
radius rvir, and took into account the uncertainties in the dis- 
tance to Draco, c/ = 82 ± 6 kpc, and its line-of-sight velocity, 
Vios = -293 ±2 km s"', which were taken from Mateo 1998.) 
One can see that the observational results are virtually incon- 
sistent with Draco moving along a bound orbit around the 
Milky Way. More quantitatively, the chance for a bound orbit 



Orbit e R 



rad 



p 
kpc 



Ra Ra/Rp 

kpc 



p 

Gyr 



V,.a 

km s"' 



fp 



0.37 70.1 
0.61 62.4 
0.85 51.1 
1.09 35.7 



260 
162 
122 
104 



1.33 18.5 96.9 
1.54 2.47 92.1 



3.7 
2.6 
2.4 
2.9 
5.2 
37 



4.86 
2.99 
2.19 
1.71 
1.42 
1.23 



78.3 
103.9 
112.1 
99.9 
65.5 
11.5 



0.25 
0.47 
0.76 
0.64 
0.57 
0.54 



Note. — Here Vi.a is the tangential velocity at R = Ra and fp is the fraction 
of the time Draco spends at the distances = 70 ... 1 40 kpc from the Galactic 
center, where 75% of Galactic dwarf spheroidals are currently located. 



with the radial period P <S Gyr is only 5.5% (or 7.4% for any 
period, which includes periods much longer than the Hubble 
time). This is not surprising, as the proper motion measure- 
ments of Scholz & Irwin ( 1994) imply that Draco moves in 
the Galactic halo with a staggering speed of 610 ± 190 km s"' 
(one sigma error bars), with no realistic Milky Way model 
being able to keep it gravitationally bound. 

Given that Draco appears to be a pretty normal dwarf 
spheroidal galaxy and that the dwarf spheroidals strongly 
concentrate toward the tw o large spi rals in the Local Group 
(Mashchen ko. Carignan. & Bouchardl l2004 . the Milky Way 
and the M31 galaxy, it seems to be very unlikely that this 
dwarf moves along an unbound orbit around our Galaxy, 
and just by chance happened to be at the present small dis- 
tance from the Galactic center. We assume instead that 
Draco moves on a bound orbit, wit h P < 8 Gyr, and that 
the proper motion measurements of 'S cholz & IrwinI il994h 
are wrong. Proper motion measurements of the outer Galac- 
tic halo objects based on heterogeneous collection of pho- 
tographic plates spanning a few decades, such as those of 
Scholz & Irwin ( 1994), are notoriously difficult to correct for 
all possible sources of systematic errors. As another ex- 
ample, the prope r motion measurement of the same authors 
JScholz & Ir wih 1994) for Ur sa Minor is more than one sigma 
away from each of the three more recent results, including 
two derived with the help of the Hubble Space Telescope 
(!Piate k et alJ2005[ their Fig. 16). 

Fortunately, even in the case when the proper motion is not 
known, some constraints can be placed on the orbital elements 
of Draco. We noticed that pericentric and apocentric distances 
for all the bound orbits from Figure |6^ are not completely 
independent, and follow closely the solution of the equation 



[Rl + (7000/RafY^^ = 76kpc, 



(22) 



which is a circle with the radius of 76 kpc in the coordinates 
X = Rp, y = 1000/ Ra (see Figure (SJi). (Here Rp and Ra are in 
kpc.) As one can see in Figure |6j3, despite the observational 
uncertainties in d and Vios and the fact that the proper motion 
in not known, the Draco's bound orbits stay in a relatively 
narrow zone near the solution of equation ( I22t . Neglecting 
the spread of the points around the circle, one can assume 
then that bound Draco's orbits form a one-dimensional family 
of models, with both Rp and Ra depending only on the polar 
angle 9, with tan0 =y/x = 1000/(RpRa). 

We chose six different values of the polar angle 9 (shown 
as radially divergent lines in Figure|6j3, and listed in Table|3 
to obtain a sequence of bound Draco's orbits. We used plots 
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Fig. 6. — (a) Proper motion vectors for Draco. The observational result of lScholz & Irwfir^i 19941) is shown as a cross with the circle representing a 0.5 mas yr"' 
one-sigma en'or bar. Dots correspond to proper motion vectors which result in a bound orbit with the period P < 8 Gyr in the Milky Way potential, (b) Draco's 
bound orbits with P < 8 Gyr in the (Pp,7000/Pa) coordinates (dots). In these coordinates, the orbits follow well a circle with a radius 76 kpc (thick solid line). 
The cutoff at 7000/Ra ~ 20 kpc is due to the imposed cutoff in the orbital period P < 8 Gyr. For six different values of the polar angle d (numbered radially 
divergent straight lines) we plot the averaged distances of the points from the reference point (circles with crosses). As one can see, all the averaged distances are 
very close to 76 kpc. 




B, rad B, rad 

Fig. 7. — Our choices for Draco's orbits. Dots show bound orbits with P < 8 Gyr. Vertical lines correspond to the six different values of the polar angle Q 
from Figure|5I). Circles with crosses mark the averaged orbital parameter (either pericentric distance Rp, panel a, or apocentiic distance Ra, panel b) values for 
different angles d. 



shown in Figure to estimate the typical values of Rp and 
Ra corresponding to the particular values of the angle Q. In 
Table 12 we list the parameters of the derived Draco's orbits. 
The range of covered orbital periods is ^ 1-5 Gyr. All the 
orbits except for the last one have a pericenter at the distances 
where the tidal disruption properties of our NFW halo are 
comparable to those of the composite Milky Way model of 
iJohnston etaHlH^ see Figure|5}. The orbit 6 was designed 
to explore the extreme case of a virtually radial orbit. The 
apocenter of the orbit with the longest period is at 260 kpc, 
which is at the very edge of the virialized Galactic halo. As 



one can see from Table |2] the derived sequence of orbits is 
not trivial, with the orbits being more eccentric for longest 
and shortest periods, and becoming rounder for intermediate 
periods of ~ 2 Gyr. 

It is interesting to note that six out of eight, or 75% of the 
"classical" Galactic dwarf spheroidals (we exclude Saggitar- 
ius as being currently tidally disrupted) are located in a rather 
narrow interval of Galactocentric distances R = 1Q... 140 kpc 
(Mateo 1998). This includes Draco, and excludes Leo I and 
II. In Table |2l we list for each orbit the fraction of time fp 
Draco spends in this interval of Galactocentric distances. Sta- 
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tistically speaking, if Draco is a "normal" dwarf spheroidal, 
fp should be around 0.75. Our orbit 3 is in this sense the like- 
liest orbit for Draco, and orbit 1 (and probably 2) are rather 
unlikely. 

In Figure |4] we show as long-dashed lines the halos with 
the analytical tidal radius rtid equal to 0.85 kpc, which is 
the radius where the den sity of Draco stars on the map of 
lOdenkirchen et all (1200 Ih is two sigma above the noise level. 
We calculated from 



'"(Aid) 



'tid 



R dM 
MiR)~dR 



M(R) 



(23) 



(|Havashi et al. 2003), where M(R) and m(r) are the enclosed 
mass for the Milky Way and the satellite, respectively. Even at 
this lar ge distance from the Draco's center Odenkirchen et aj 
^OT) did not see any sign of Draco being tidally distorted by 
the Milky Way tidal field. From Figure 0]one can see that all 
our best-fitting models (solid circles) have tidal radii larger 
than 0.85 kpc (even for the worst orbit with Rp = 2.5 kpc). 
One would naively conclude that the tidal forces are not im- 
portant for our Draco models. But the reality is more compli- 
cated than that. Numerical A^-body simulations of the evolu- 
tion of subhalos orbiting in the host halo showed that removal 
of DM from the outskirts of the satellite results in the expan- 
sion of the satellite, which reduces its average density and ex- 
poses more DM to the action of the tidal field dHavashi et al.l 
12003; Kazantzidis et al. 2004). Burkert halos have lower av- 
erage density than NFW halos, and as a result a re even eas- 
ier to disrupt tidally iMashchenko & Sillsll2005bh . Moreover, 
even well inside the tidal radius, distribution of bound stars 
can be noticeably distorted by the tidal filed of the Milky 
Way, which would be at odds with the Draco observations of 
iQdenkirchen et al. (2001 ). To correctly describe the above ef- 
fects, we had to resort to high-resolution A^-body simulations 
of the tidal disruption of our best-fitting composite (DM + 
stars) models from Table^in the static potential of the Milky 
Way. This will be described in the following sections. 

4.2. Isolated Models 

To run the A^-body simulations of the tidal strip- 
ping of Draco in the static gravitational potential of the 
Milky Wa y, we used the parallel tree-code Gadget- 1.1 
dSpringel. Y oshida & White 2001 ). We generated equilibrium 
DM halos for our eight models (see Tabled using the pre- 
scription in Mashchenko & Sills (2005a). The essence of this 
approach is to use explicitly the distribution function (DF) to 
set up the initial distribution of velocity vectors of DM parti- 
cles (' isotropy was assumed). It was argued ( Kazantzid is et al.l 
120041) that using DFs explicitly is far superior to traditional 
local Maxwellian approximation for cuspy models such as 
NFW. To reduce the boundary effects, we truncate DM ha- 
los at a distance of two virial radii rvii from the center This 
results in virtually no evolution for isolated models within 
the virial extent of the halos after 10 Gyr We chose the 
gravitational softening lengths (separately for DM and stars) 
to be commensurable with the average interparticle distance: 
e = 0.77r/,A^"'/^ ( Havashi et alJ2003h . Here r/, is the half-mass 
radius of the system. For stellar particles, = 8 .6 pc. For DM 
particles, Edm values are listed in Tabled 

To set up the initial distribution of stars inside the DM halo, 
we use the same pseudo-Maxwellian approximation we used 
to measure the projected line-of-sight velocity dispersion pro- 
files in §|3l(step 3). The only difference is that now we use a 



larger number of stellar particles, A^, = 30,000, to reduce the 
Poisson noise in the observable properties of the models. 

Our baseline number of DM particles is A^dm = lO*". Test 
runs showed that in the most massive models (Nl,2 and B 1,2) 
a much larger number of DM particles is required to prevent 
the artificial evolution of the stellar cluster at the center of 
the halo. We observed the central stellar density being signif- 
icantly reduced (by more than an order of magnitude in the 
worst cases) in the low-resolution models. This effect does 
not depend on the number of stellar particles, and becomes 
less severe for larger A^dm and/or £dm- As stars are only a 
trace population in our models, the most obvious explana- 
tion for the above artifact is that stars get scattered from the 
imperfections of the granulated gravitational potential of the 
DM halo. Setting A^dm = lO*" for the models Nl and Bl, and 
Ndm = 3.16 X 10^ for the models N2 and B2 resulted in ac- 
ceptable level of the artificial evolution of the surface bright- 
ness profile I](r) of the stellar cluster. In all our isolated mod- 
els, the change in the central surface brightness was negli- 
gible after 10 Gyr of evolution (with the only exception of 
the model Bl, where the change was -0.4 dex). The radius 
corres ponding to the outmost r eliable Draco isodensity con- 
tour of lOdenkirchen et alJ (12001 ) with E = 12,060 Mq kpc'^ 
(r = 0.85 kpc initially) increased by mere 0.03 - 0.05 dex. 

Unfortunately, we observed significant evolution in the ve- 
locity anisotropy profiles ri(r) in our isolated models, espe- 
cially at the center of the halo. All our best-fitting stellar mod- 
els have a strong radial anisotropy at the center (see Table^. 
In our A^-body simulations with live DM halos the stellar core 
becomes close to isotropic within the first Gyr, suggesting that 
the particle-particle interactions are not the main culprit, as 
such effects would take Gyrs to manifest themselves. This 
became more apparent after we ran simulations for all our 
models from Tabled with A^* = 30,000 stellar particles and 
static DM potential. In this setup, close encounters between 
particles are extremely rare due to very low stellar density. In 
the static models we see the same effect (of slightly smaller 
magnitude) as in the runs with live DM halos: central radial 
anisotropy reduces almost to zero after a few crossing times. 
We do not know the exact reason for the above effect. The 
possible explanations are (1) the (unknown) DF correspond- 
ing to our choice of anisotropy, stellar density, and gravi- 
tational potential profiles is okay (positive everywhere), but 
the pseudo-Maxwellian approximation we use to set up the 
stellar velocities is not good for systems with strong radial 
anisotropy at the center, and/or (2) the DF is not physical (i.e. 
is negative) at the center. We also want to point out that a real 
stellar system cannot have a radial anisotropy all the way to its 
center, as the radial velocity dispersion diverges in such case 
when r — > 0. To demonstrate this, we write down the solution 
of the Jeans equation il5l in the case of constant anisotropy 
(lLokasl2001l) 



1 



dr 



(24) 



where the constant C, = 2(3 = Arj / { \ + rf) is positive for the case 
of radial anisotropy. The integral in equation \2Al is non- 
zero for r = 0, resulting in divergent cr,- at the center of the 
system. Obviously, in a real object radial anisotropy should 
break down at some radius, changing into isotropy or tangen- 
tial anisotropy. It is also hard to expect radially divergent ve- 
locity dispersion in a stellar system with a flat core, especially 
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Fig. 8. — Line-of-sight stellar velocity dispersion profiles for Draco for the 
models evolved in a static DM halo potential. The notations are the same as 
in Figure|5] 



in the case of a flat-cored (Burkert) DM halo. 

We want to emphasize that even though we cannot guar- 
antee that the stellar models in Table [2 are physical (espe- 
cially at the center), in our simulations they quickly relax 
to a stable configuration, with the surface brightness profile 
virtually identical to the observed one, and the line-of-sight 
velocity dispersion profile trios still consistent with the ob- 
servations. In Figure|8lwe plot crios('") profiles for our stellar 
models in a static DM potential after 10 Gyr of evolution inte- 
grated over the same six radial bins as the observational data 
of Wilkinson et al. (2004). As one can see, the profiles for 
evolved stellar clusters are a reasonably good match to the 
observations. (The x' measure becomes a factor of two larger 
in the evolved models.) 

We also used the isolated runs to estimate the accumulated 
total energy errors in our models. In the case of live DM halos, 
the total energy for DM + stars is conserved to better than 
--0.1% (typically - 0.05%) after 10 Gyr of evolution. To 
estimate the total energy errors for stars only, we measured the 
difference in stellar total energy at f ~ 3 Gyr (when the cluster 
has reached a steady state configuration) and at f = 10 Gyr 
in our static DM halo simulations. The difference was again 
<0.1%. 

4.3. Tidal Stripping Simulations 

We simulated evolution of the eight stars+DM models of 
Draco from Table orbiting for 10 Gyr in the static spheri- 
cally symmetric potential of the Milky Way given by equa- 
tion ( I19> . To run the simulations, we us ed the parallel versio n 
of the multi-stepping tree code Gadget (Springel et al. i200ll) . 
The number of stellar and DM particles and the corresponding 
gravitational softening lengths were the same as in the isolated 
models described in § 14.21 Each Draco model was simulated 
for six different orbits from Table |2l Initially, Draco was lo- 
cated at the apocenter of its orbit. Altogether we made 48 
tidal stripping simulations. We refer to the models by both 
the model number from Table ^ and the orbit number from 
Table|2j e.g. model Nl-5. Most of the results we give below 
are for the latest moment of time when the dwarf was located 



at the current distance of Draco from the Galactic center of 
~ 82 kpc (we do not make distinction between the dwarf mov- 
ing inward or outward), which took place 7.4-10 Gyr from 
the beginning of the simulations (9.4-10 Gyr for the orbits 
3-6). 

All our models experienced tidal stripping of different de- 
gree. As Figure shows, the DM mass of the gravitation- 
ally bound remnant is between 90% (model N5-1) and 0.1% 
(model Nl-6) of the original mass at the end of the simula- 
tions. Three of our models became completely gravitationally 
unbound within 10 Gyr: Bl-6 (after 3.2 Gyr and 3 pericenter 
passages), B2-6 (after 7.8 Gyr and 6 pericenter passages), and 
Bl-5 (after 9.4 Gyr and 7 pericenter passages). In the latter 
case, the dwarf is still bound when it is located at ~ 82 kpc 
at the end of the simulations (when we compare the results 
of the simulations with the observed properties of Draco) - 
but only barely so. It is not surprising that all the unbound 
models have Burkert halos. Indeed, for given mass and scal- 
ing radius these halos have lower averaged density than NEW 
halos in the central area. If truncated instantaneously at a cer- 
tain radius, the total energy of the remnant becomes positive 
for Bu rkert halos at a radius ^2.1 times larger than for NEW 
halos ( Mashchenko & SiUsl2005"bl) . 

Another result which can be explained is that the most mas- 
sive halos are easier to strip and disrupt tidally than the less 
massive ones. All our models (both massive and of lower 
mass) have a comparable DM density within the observed 
extent of Draco (because of the virial theorem), so in the 
point mass approximation they should be equally susceptible 
to tidal forces. However, the point mass approximation (used 
to derive equation 1231 ') breaks down for our most massive 
halos, as their size becomes comparable to the pericentric dis- 
tance, so the strongly non-linear components of the tidal force 
become important. Not surprisingly, our most massive halos 
on the orbit 6 were either totally disrupted (model B 1 -6) or 
lost 99.9% of its original mass (model Nl-6) after 10 Gyr of 
evolution. 

In all of our models a fraction of stars has been tidally 
stripped by the end of the simulations - even for the orbit 1 . In 
two cases (models Bl-6 and B2-6), stars have become com- 
pletely unbound by the time the dwarf was passing at the dis- 
tance of ^ 82 kpc from the Galactic center for the last time. In 
other cases, the fraction of escaped stars was between ^ 10"^ 
for the models B 1-2,3,4 and 96.6% for the model Nl-6 (see 
Eigure|9j)). 

When analyzing the global properties of the stellar cluster 
at the end of the simulations, the most obvious result is the 
fact that the more disruptive the orbit is, the more the system 
is affected by the tidal shocks experienced near the pericen- 
ter of the orbit. (Orbits with a larger number from Table |2l 
are more disruptive for two reasons: they have smaller peri- 
centric distance, and they have shorter orbital period, so the 
number of pericentric passages in 10 Gyr is larger.) Dwarfs 
are puffed up by the tidal shocks, with the central line-of-sight 
velocity dispersion (Tq becoming smaller (Eigure|9j;), central 
surface brightness decreasing (not shown, but is qualitatively 
similar to (Tq behavior), and the projected half-light radius be- 
coming larger (except for the extreme case of the orbit 6 when 
both tidal shocking and tidal stripping are important, see Fig- 
ure|9ll). 

Our models are highly idealized when it comes to long-term 
evolution of stellar tidal streams. On our most disruptive or- 
bits 5 and 6, a significant fraction of stars becomes unbound 
over the course of the tidal evolution, with the most of the 
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Fig. 9. — Draco models' parameters near the end of the simulations (when the dwarf was ~ 82 kpc away from the Galactic center) as a function of orbit. 
Sohd/dashed hnes correspond to models with NFW/Burkert DM halo profiles, respectively. Thick lines con'espond to the most massive halos (models Nl and 
Bl). (a) Gravitationally bound DM mass Mdm- (b) Fraction of stars having become unbound /, . (c) Central line-of-sight stellar velocity dispersion (jq. (d) 
Projected half-light radius for the bound stellar cluster r;, . To facilitate the comparison of different models, both ctq and r/, are normalized to the same value for 
the orbit (corresponding to isolated models). 



tidally stripped stars following very closely the almost radial 
orbit of the dwarf. From the Sun location, many of these stars 
project on a rather small area in the sky inside or around the 
apparent location of the dwarf. This can be quite unphysi- 
cal, as such cold tidal streams are not expected to survive for 
many gigayears in the Milky Way halo due to its triaxiality 
and dumpiness (as predicted by cosmology) and due to in- 
teraction with baryonic structures in the Galaxy (stellar disk 
with spiral arms, stellar bar, giant molecular clouds). To cir- 
cumvent this difficulty, in this section we discuss all the ob- 
servable model properties for two extreme cases: (a) all stars 
are taken into account, and (b) only the stars (both bound and 
unbound) located within the spatial distance of 5 kpc from the 
center of the dwarf are considered. In the latter case, the spa- 
tial truncation of the tidal stream ensures that only the most 
recently stripped stars are used for calculating the observable 



properties of the models. 

On a more detailed level, the impact of both tidal shocks 
and tidal stripping and heating can be seen in Figure^] Here 
we show the surface brightness profiles (panel a) and line- 
of-sight velocity dispersion profiles (panel b) for the models 
Nl-1,4,5,6. The most obvious effect in FigureFTOlis the global 
decrease of ctio.s for orbits with smaller pericentric distance 
(excluding the orbit 6), accompanied by a small decrease in 
the central surface brightness and a slight radial expansion of 
the system. No obvious tidal truncation and tidal heating is 
observed in the outskirts in the cluster. These results were 
obtained for the stars located within the spatial distance of 
5 kpc from the center of the dwarf, but in the case when all 
stars are included the profiles are practically the same. The 
orbit 6 is a completely different case: one can see the trios 
being dramatically inflated in the outskirts of the dwarf due 
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Fig. 10. — Observable properties of the model Nl for a few different orbits near the end of the simulations (when the dwarf was ^ 82 kpc away from the 
Galactic center). Only stars (both bound and unbound) located within the spatial distance of 5 kpc from the dwarf's center are considered. SoHd, dotted, short- 
dashed, and long-dashed lines correspond to orbits 1, 4, 5, 6, respectively. T he obs erver is assumed to be located at the Galactic center, (a) Surface brightness 
profile S(r). We also show the observed Draco profile of Odenkirchen et aT] l200H their sample S2). (b) Stellar line-of-sight velocity dispersion profile cr\a^(r). 
We also show the observed Draco velocities o fi Wilkinson et aU ii2004V 



to superposition of tidally removed stars on the dwarf (Fig- 
ure ^^). This effect becomes even more dramatic when we 
include all the stars, with the line-of-sight velocity dispersion 
reaching 50-70 km s"' in the outermost observed bin. As 
we discussed in the previous paragraph, it is quite unlikely 
that old tidal streams can stay so well collimated for many 
gigayears to produce the above effect. But even for the con- 
servative case of considering only freshly stripped stars the 
steeply growing dies profile for the model N 1-6 is grossly in- 
consistent with the observed profiles of Draco and other dwarf 
spheroidals, where (Tio,s is either not changing or decreasing at 
large distances from the center 

The change in the surface brightness profile for the model 
Nl-6 is also quite dramatic (Figure [Toh. long-dashed line). 
The outer S(r) slope becomes very shallow which is inconsis- 
tent with the observed profiles for Galactic dSphs. The slope 
becomes even more shallow when we include stars beyond the 
spatial distance of 5 kpc from the center of the dwarf. Similar 
behavior (in terms of shallow outer S profiles and inflated (Tios 
in the outskirts of the dwarf) is also observed in other models 
with the orbit 6 (both bound and unbound). The possible ex- 
ceptions are the models N4-6 and N5-6, which are consistent 
with the observations of Draco when we consider only freshly 
stripped stars. We conclude that it is unlikely that Draco and 
other Galactic dSphs have experienced tidal interactions as 
dramatic as our models on the orbit 6. 

To facilitate the comparison of our models with 
the obser ved ste llar isodensity contours of Draco of 
lOdenkirch en et alJ (|2001) we designed the following pro- 
jection algorithm, (a) The frame of reference is rotated to 
place the center of the dwarf on the negative side of the axis 
Z, with the axis X located in the orbital plane of the dwarf 
and pointing in the direction of the orbital motion, (b) As 
the direction of the proper motion of Draco is not known 
(see discussion in § 14. U and the current angle between the 
vector connecting Draco with the center of the Galaxy and 
the vector connecting Draco with Sun is c/? = 5.°7, we consider 



three extreme cases of the possible vantage point location 
which should encompass the whole range of projected model 
appearances: view from the Galactic center (jp = 0°), view 
from a point in the orbital plane of Draco located at 82 kpc 
from the dwarf with = 5.°7, and view from a point in the 
plane which is perpendicular to the orbital plane of Draco 
at 82 kpc and with Lp = 5?7. We found that due to the fact 
that for Draco the angle ip is small (which is also the case 
for other Galactic dSphs with the exception of Sagittarius), 
the observable properties of our models (S and crios profiles 
and surface brightness maps) do not depend noticeably on 
the vantage point we choose - especially for the case when 
we only include freshly stripped stars, (c) We exclude stars 
with z > -8.5 kpc to avoid contamination of our maps with 
local tidal stream overdensities which would be discarded by 
observers as local Galactic stars, (d) We perform prospective 
projection of the stellar particles smoothed with a Gaussian 
beam which has a fixed physical size (either 0.15 or 0.5 kpc) 
and brightness inversely proportional to the square of the 
distance of the particle from the observer This procedure 
makes the "surface brightness" of individual particles invari- 
ant of the distance from the observer, which is appropriate for 
spatially resolved clumps of the tidal stream. 

The most interesting result obtained from the analysis of the 
surface brightness maps is the lack of the classical S-shaped 
tidal tails in our models. In the case of significant tidal strip- 
ping (Figure ^2 two different vantage points are shown) the 
stellar isodensity contours change from being spherical near 
the center of the dwarf to being increasingly more elliptical 
and often off-centered at larger distances. In the cases with 
less severe stripping, outer contour boxiness is observed in 
some of the models. For orbits 1-3 the surface density of 
the tidally removed stars is so low that it is hard to draw 
any conclusion as to the shape of the tidally distorted isoden- 
sity contours (except for the fact that the galaxy is observed 
against the background and/or foreground of a few degrees 
wide belt of tidally stripped stars). The explanation for the 
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Fig. 1 1 . — Stellar surface brightness maps for the model N2-6 near the end of the simulations (when the dwarf was ~ 82 kpc away from the Galactic center). 
The contour levels are (8, 16,32,... ,512) X 10^ M0 deg^ (calculated for the nominal distance of Draco of 82 kpc). Cross m arks the center of the dwaii'. Thin 
line circle has a radius of 0.85 kpc - the radius of the Draco two-sigma surface brightness contour of lOdenkirchen etall 120011) . In the upper left corner we show 
the size of the Gaussian beam used to make the maps, (a) The observer is in the plane of the Draco's orbit, (b) The observer is in the plane perpendicular to the 
Draco's orbit. 



above effect is in the facts that the tidal stripping is signifi- 
cant only for very eccentric (almost radial) Draco orbits with 
Rp < 20 kpc, and that currently Draco is located > 10 kpc 
away from the apocenter of its orbit (see Table |3- Under 
these circumstances, line of sight practically coincides with 
the direction of the tidal elongation of Draco, with both tidal 
tails seen edge-on. This is an interesting result, as it suggests 
that the lack of S-shaped isodensity contours in the outskirts 
of the Galactic dSphs cannot be used to support a claim that 
the dwarf has not experienced significant tidal stripping in the 
gravitational potential of the Milky Way. It appears that the 
presence (absence) of tidally heated stars in the outskirts of 
dSphs is much more sensitive indicator of tidal stripping be- 
ing significant (not significant). 

We measured for all our models the critical surface bright- 
ness Ec when the outer isodensity contours become notice- 
ably distorted due to tidal effects. The two models which 
became gravitationally unbound by the end of the simula- 
tions (Bl-6 and B2-6) show very distorted contours all the 
way to the center of the dwarf (see Figure [T2I . They also 
have dramatically inflated line-of-sight stellar velocity dis- 
persion profiles (up to 50-100 km s"'), and very shallow 
outer surface brightness slopes. All of the above makes these 
two models (and most probably any other unbound model for 
a dSph) completely inconsistent with the observed proper- 
ties of the Galactic dSphs. Among the bound models, the 
one which has experienced the most dramatic tidal stripping 
(model Nl-6) is also the one which has the largest value 
for 'Sc- ~ 4 X lO'* M0 deg"^, which would correspond to 
~ 3 sigma isodensity contour of lOdenkirchen et"al] J2001I) . 
In all other cases, the value of is significantly lower: 
(7.4- 17) X 10^ M0 deg"^ for the models on the orbit 6, and 
less than 10^ Mq deg"^ for other orbits. For the orbits 1, 2, 
and 3 E^. was too small to be measured. These values of E^ are 
significantly lower (by a facto r of > 1.5) th an the two-sigma 
detection limit of Odenkirche n et alJ il200l . The angular ra- 
dius of the smallest distorted isodensity contour is ~ 1° for 
the orbit 6 (with the exception of the models Bl, B2, and Nl), 
and > 1 .°8 for other orbits. 

One very interesting special case is that of the model Bl. 



This model has a halo with a flat DM core of size r, ^ 1 .4 kpc 
(see Tabled), so all the observed extent of Draco is within this 
large harmonic core. In Figure we show the inner (corre- 
sponding to the observed part of Draco) isodensity contours 
for the models Bl-4 and Bl-1. Unlike all other models (both 
NFW and Burkert), here one can see a relatively strong tidal 
distortion of the contours at distances < 0.°5 from the center 
of the dwarf. The distortions are substantial even for orbit 
1 (Figure 113b ) which has the largest pericentric distance of 
70 kpc. The distorted contours are both elliptical and non- 
concentric. Interestingly, the effect is minimal for the two 
orbits which have the smallest eccentricity - orbits 2 and 3 
with Ra/Rp ^ 2.6. The distortions are strong for the more ec- 
centric orbits 1, 4, and 5. It appears that it is the variability 
of the tidal force rather than its strength which is the primary 
governing factor for the above effect. As the observed stellar 
isodensity contours of Draco are very regular and concentric 

( Odenkh'chen et al 2001 ) . results of our simulations suggest 

that it is unlikely that Draco resides in a large DM harmonic 
core - unless it happened to move on a close to circular orbit 
with Rp ~ 45 -65 kpc and Ra/Rp< 2.6. 

5. DISCUSSION 

In this paper we presented a sequence of composite (stars 
+ DM) models for Draco, listed in Table [0 which satisfy 
all the available observational and cosmological constraints. 
We showed that for the most of physically plausible orbits 
of Draco in the Galactic potential the tidal forces could not 
modify the observable properties of our models appreciably 
after 10 Gyr of evolution. Both "standard" cuspy NFW DM 
halos and Burkert halos with a flat core provide a reasonably 
good description of Draco. The properties of a Burkert halo 
are better constrained by our analysis. Most importantly, if 
Draco has a flat core, it should have formed at or before the 
end of the reionization of the universe: z > 6.5. Tidal strip- 
ping simulations put even stronger constraints on the flat-core 
case: we showed in the previous section that our most massive 
Burkert model Bl would be consistent with the observations 
only for a very limited range of possible Draco orbits: orbit 6 
is ruled out as the model becomes completely unbound with 
dramatically inflated crios and very shallow E profiles, whereas 
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Fig. 12. — Surface brightness maps for the two unbound models at the end of the simulations (when the dwaif is located at /? ~ 82 kpc). Only freshly stripped 
stars (within a radius of 5 kpc from the dwarf) are included. The observer is located at the center of the Galaxy. Cross marks the center of the dwarf. In the upper 
left comer we show the size of the Gaussian beam used to make the maps, (a) Model Bl-6. Contours are (4. . . 19) X 10^ Mq deg"^. (b) Model B2-6. Contours 
are (1...8.6) X lO' Mq deg-^. 




for the orbits 1, 4, and 5 (and also 6) we observe significant 
distortion of inner stellar isodensity contours which is most 
definitely not consistent wi th the regular isodensity c ontour 
shape observed in Draco bv lOdenkirchen et alJ J2001h . Only 
the lowest eccentricity orbits 2 and 3 are not ruled out by our 
analysis. 

An NFW halo is less constrained by the available observa- 
tions of Draco: the halo formation redshift z is anywhere be- 
tween ^ 2 and ^10, whereas the initial virial mass could be 
between ^10^ and 5 x 10^ Mq. In the smaller mass (and 
larger z) limit the halos are so sturdy that even for our most 
disruptive orbit 6 the observable parameters of the model can 
still be consistent with the Draco observations after 10 Gyr of 
tidal evolution. If Draco was accreted by the MiUcy Way more 
recently than 10 Gyr ago, the impact of the tidal forces would 
be even smaller For more massive NFW halos and for all 
Burkert models the orbit 6 can be ruled out as the model pre- 
dicts inflated line-of-sight velocity dispersion in the outskirts 
of the dwarf which would be at odds with observations. 



How strong is our case against Draco (and other dSphs) be- 
ing a "tidal dwarf" - remnants of a dwarf galaxy which are 
not gravitationally bound at the present time? In the previ- 
ous section we suggested that the fact that the line-of-sight 
velocity dispersion is dramatically (by up to an order of mag- 
nitude) inflated in our two "tidal dwarf" candidates, models 
Bl-6 and B2-6, can be used to rule out the "tidal dwarf" ex- 
planation for Draco. Here we want to caution that a more 
detailed comparison between the model and observations is 
required to critically assess our conclusion. Our large esti- 
mates of (Tios were derived for stars with any line-of-sight ve- 
locity projected onto the dwarf disk and optionally restricted 
to lie within the spatial distance of 5 kpc from the dwarf's 
center. Many of the high-velocity tidal tail stars responsible 
for inflating cios would be discarded by observers as "not be- 
longing to the galaxy". In FigureOwe show the distribution 
of line-of-sight velocities Vios for the models Bl-6 and B2- 
6. We show separately histograms for freshly stripped stars 
(solid lines) and for all stars projected on the disk of the dwarf 
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Fig. 14. — Distribution of line-of-sight velocities Vios for the models which became unbound by the end of the simulations. The observer is located at the center 
of the Galaxy. The dwaif is located ~ 82 kpc away from the center of the Galaxy. SoHd and dashed lines correspond to the stars within the spatial distance of 5 
and 30 kpc from the densest part of the unbound dwarf, respectively. Only stars located within 0? 5 from the densest part of the galaxy are taken into account, (a) 
Model Bl-6. (b) Model B2-6. 



(dashed lines). As one can see, the situation depends strongly 
on how recently the dwarf became unbound, and on longevity 
of the cold tidal streams in the Galactic halo. The model Bl-6 
became unbound many orbits (almost 7 Gyr) ago, and has a 
very wide distribution of Vios - even for freshly stripped stars. 
The model B2-6, on the other hand, became unbound only 
~ 2 Gyr ago, and has more complex distribution of Vios- In 
this model, the freshly removed stars are virtually all concen- 
trated within a relatively narrow interval, with ctios being in- 
flated mainly due to the presence of one high velocity stellar 
particle with Vios — 160 km s"'. Such stars will definitely be 
discarded by observers. When we consider all stars (dashed 
line in Figure[T4b). the distribution of Vios is much wider than 
in Galactic dSphs. In the case of the model B2-6, the situation 
thus sensitively depends on how long the tidal stream can stay 
collimated in the Milky Way potential. 

Another important evidence against Draco being an un- 
bound stream of stars is presented in Figure [21 Here we 
show surface brightness maps for our two unbound models 
Bl-6 and B2-6. One can see that the contours are irregular 
and not concentric - even near the center of the dwarf. This is 
in sharp contrast with the regular appearance of the observed 
isodensity contours in Draco (Odenkirchen et al. 2001). 

Pending the arrival of accurate proper motion measure- 
ments for Draco, let us be slightly more definitive in trying to 
determine the nature and cosmological significance of Draco 
by assuming that it is moving on the orbit 3, which is the most 
probable one (see § 14. 1> . From Figure |9] one can then infer 
that if Draco is a cosmological halo, its current DM mass is 
between 7x10^ and 3 x lO** Mq, the fraction of the tidally 
stripped stars is < 3%, and the central line-of-sight velocity 
dispersion ao, central surface brightness Eq, and the half-light 
radius r/, has changed due to tidal shocks by no more than 
-0.07, -0.04, and 0.03 dex, respectively, in the last 10 Gyr. 
This orbit has Rp = 51.1 kpc, placing it well outside of the 
Galactic disk. Stellar tidal tails produced by our models on 
this orbit are extremely weak, with the surface brightness sen- 



sitivity required to see the isodensity contours distorted due 
to tidal forces being better than ^ 200 Mq deg~^, or more 
that two orders of magnitude better that the observations of 
'Oden kirchen et a l. (lOOl ). The DM halo could be either NFW 
or Burkert, and was formed after z ~ 1 1 . Our results then 
support either of the two re cently proposed solutions to the 
"miss ing satellites" problem (iKlvpin et al.ll999tlMoore et 
Il999h : that the Galactic dSphs are the most massive subhalos 
predicted by cosmological simulations to orbit in the halo of 
a MiU cy Way sized galaxy (Stoehr et al. 2002; Havashi et ^ 
120031) . or that the Galactic dSphs are the halos which man- 
aged to form stars before the reionization of the universe was 
completed around z ^ 6.5 (Bullock et al] l2000h . Our analy- 
sis suggests that unfortunately there is not enough of obser- 
vational data yet to discriminate between the two above sce- 
narios. Much better quality line-of-sight velocity dispersion 
profiles, deeper surface brightness maps, and ideally accurate 
proper motion measurements are required to produce further 
progress in this direction. 

We would like to mention a few most important deficiencies 
of our tidal stripping model. The first one is due to our use of 
a spherically symmetric potential for the Milky Way. As a 
result, we ignore disk shocking, which can be very important 
for the orbits with Rp < 20 kpc. We tried to partly circum- 
vent this deficiency by considering an orbit with extremely 
small pericentric distance: orbit 6 with Rp ~ 2.5 kpc. Ideally, 
we would prefer to include the disk in our simulations, but 
this would result in significant increase in number of required 
models, which would make our project not feasible with the 
current level of computing power. 

The second problem is generic to existing tidal stripping 
simulations of Galactic satellites (dSphs and globular clus- 
ters), and is caused by our use of a static potential for the 
Milky Way halo. In a realistic (live) Galactic halo very mas- 
sive satellites should experience dynamical friction, which 
would gradually bring them closer to the center of the Galaxy. 
Unfortunately, we could not use a simple analytical formula 
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to estimate the impact of the dynamical friction on our results. 
The dynamical friction would be strongest for the orbits 6 and 
5, which have the smallest pericentric distances. Our subhalos 
on these orbits experience dramatic mass loss (up to 90% and 
70%, respectively) during the first pericentric passage, render- 
ing the constant satellite mass formula of Binney & Tremainl 
Jl987h not applicable. The dyna mical friction equation of 
IColpi. Maver. & Governatol Jl999h does not have this limita- 
tion, but it was derived for the special case of a subhalo with a 
truncated isothermal DM density profile with a core which is 
very different from both the NFW and Burkert profiles of our 
subhalos. We want to emphasize that an inclusion of dynami- 
cal friction in our model would make our conclusion, that the 
observable properties of the most of our dwarfs were not af- 
fected noticeable by tidal forces, even stronger, as the dwarfs 
would start off at larger distance from the Milky Way center 
where the tidal forces are weaker. 

In addition, the use of static potential ignores the impact of 
the gravitational field of the satellite on the Milky Way halo. 
This effect can be very important for massive dwarfs on al- 
most radial orbits, with the potential of both the satellite and 
the Milky Way center violently fluctuating during the peri- 
centric passage, leading to an exchange of energy between 
the dwarf and the Galactic halo (similar to the mechanism of 
violent relaxation). The above effect would probably be im- 
portant only for the orbit 6, which we were able to rule out for 
the most of our Draco models. 

It is important to mention that our models do not cover all 
possible initial configurations of Draco. In a more general 
case, one would have to start with arbitrary initial stellar and 
DM density profiles (with the initial stellar velocity disper- 
sion profile following from eq. |6l). After 10 Gyr of evolu- 
tion in the Galactic tidal field both profiles could become sub- 
stantially different, with the line-of-sight velocity dispersion 
either increasing (due to the projection of unbound stars) or 
decreasing (due to tidal shocks; see Fig.llO>. The more gen- 
eral case would require a dramatic increase in supercomputing 
time, which would make our approach infeasible. We want 
to emphasize that despite the fact that our models probably 
do not include all possible initial Draco configurations, they 



do constitute a family of fully self-consistent models which 
match well all the available observations of Draco. 

A potentially important evolutionary factor not included in 
our model is an interaction between Draco and dark subha- 
los, predicted to be present in the Milky Way halo in large 
numbers by ACDM cosmological mode ls. Thi s effect was 
studied on larger scales by Moore. Lake. & Katzl (fl998). who 
showed that in a cluster environment the galaxy-galaxy ha- 
rassment can be substantial. It is not clear if the harassment 
on a smaller, galactic scale would be of the same order: unlike 
the cluster scale, where the numbers of modeled and observed 
galaxies are in good agreement, there is a "missing satellites" 
issue on galactic scales. 

6. CONCLUSIONS 

We presented two one-parameter families (separately for 
NFW and Burkert DM density profiles) of composite (stars 
+ DM) models for Draco which satisfy all the available obser- 
vational and cosmological constraints. We showed that these 
models can survive tidal shocks and stripping on most real- 
istic Draco orbits in the Galactic potential for 10 Gyr, with 
no appreciable impact on their observable properties. Both 
NFW and Burkert DM halo profiles are found to be equally 
plausible for Draco. The DM halos are either massive (up to 
~ 5 X 10^ Mq) and recently formed (z ~ 2 ... 7), or less mas- 
sive (down to ~ 7 X 10^ M0) and older (z ~ 7 ... 1 1). Con- 
sequently, our results can be used to support either of the two 
popular solutions of the missing satellites problem - "very 
massive dwarfs" and "very old dwarfs". Higher quality ob- 
servations (line-of-sight velocity dispersion profiles, surface 
brightness maps, proper motion measurements) are required 
to further constraint the properties of Galactic dSphs and to 
place them in the right cosmological context. 

We would like to thank Jan Kleyna for providing the ob- 
served line-of-sight velocity dispersion profile for Draco. The 
simulations reported in this paper were carried out on McKen- 
zie cluster at the Canadian Institute for Theoretical Astro- 
physics. 



APPENDIX 

DERIVATION OF SPACE VELOCITY VECTOR COMPONENTS FOR GALACTIC SATELLITES 

In this section we derive the components of the space velocity vector of an object with known distance from the Sun D, 
heliocentric line-of-sight velocity Vios, and two proper motion components jia cos 5 and /i^. We work with the frame of reference 
where the center is at the Sun, the axis X is directed toward the Galactic center, the axis Y is pointing at (/ = 90°, b = 0), and 
the axis Z is directed toward the North Galactic Pole {b = 90°). Here {l,b) are the galactic coordinates. The frame of reference 
is at rest rel ative to the Galactic center. The Solar velocity vector in this frame of reference is V© = {10.0,225.25,7.17} km s"' 
dPehnen & Binnev 1998) . We assume that the Sun is located at the distance Rq = 8.5 kpc from the Galactic center. 

To convert the proper motion vector components from equatorial to galactic frame of reference, one can use 

^/cosZ7 = (^QCOS(5)cos(/?-^<5sini^, (Al) 

Mi = (Ma cos ^) sin + cos (/5 , (A2) 

where the angle tp is obtained from 

cos iy9 = (sin (5ngp - sin 5 sin b) / (cos 5 cos b) , (A3) 

sin 1^ = - sin(a - aNCp) cos 5ngp / cos b. (A4) 

Here {a^S) and {l,b) are the equatorial and galactic coordinates of the object, respectively, and (aNGP,<^NGp) are the equatorial 
coordinates of the North Galactic Pole (for the 72000 equinox, qngp = 192.°859 and 6ncv = 27.° 128). 

One can show that the three components of the space velocity vector of the object in the frame of reference moving with the 
Sun are 
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U, = zVios+Di^hix'+f)'^\ (A5) 
= [xiVios-zU,) -D(/i, cos b)yix^ +y^)'^^]/ix'+y^), (A6) 
Uy = [Difx, cosbXx^ +/) 1/2 +yu,yx, (A7) 

where x = cos I cos b, y = sin Z cos/?, and z=smb are the components of a unit vector directed from the Sun toward the object, and 
the units for D, U, and the two proper motion components {^icosb , jjih) are km, km s"\ and rad s"', respectively. 

In the frame of reference which is at rest relative to the Galactic center, the space velocity vector of the object is V = U+ Vq. 
In the cylindrical Galactic frame of reference, the three components of the space velocity vector of the object are 

Il = {S,V,+SyVy)/{Sl+Slfl^, Q = {SyV,-S,Vy)/{Sl+Slfl\ W = V„ (A8) 

where 11 is directed outward from the Galactic center in the plane of the Galaxy, 8 is the circular rotation speed in the plane of the 
Galaxy (positive for the Sun), W is directed toward the North Galactic Pole, and S = {xD-RQ,yD,zD} is the vector connecting 
the Galactic center with the object. In the spherical Galactic frame of reference, the radial and tangential velocities of the object 
are 

K = (V-S)/|S|, y, = (|v|2-y2)V2. (A9) 
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